Glassy Dynamics of Protein Folding 
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A coarse grained model of a random polypeptide chain, with only discrete torsional degrees of 
freedom and Hookean springs connecting pairs of hydrophobic residues is shown to display stretched 
exponential relaxation under Metropolis dynamics at low temperatures with the exponent /3 ~ 1/4, 
in agreement with the best experimental results. The time dependent correlation functions for 
fluctuations about the native state, computed in the Gaussian approximation for real proteins, have 
also been found to have the same functional form. Our results indicate that the energy landscape 
exhibits universal features over a very large range of energies and is relatively independent of the 
specific dynamics. 
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A huge amount of effort has recently been invested 
in modeling the interactions responsible for yielding the 
native states of proteins as their thermodynamic equi- 
librium state It has recently begun to be appre- 
ciated that such features of real proteins as the density 
of vibrational energy states || may be reproduced by 
coarse-grained model hamiltonians which capture the es- 
sential mechanism driving the folding process, namely 
hydrophobic interactions |5]-|)j. In this paper we intro- 
duce and study a model of N coupled, over-damped tor- 
sional degrees of freedom with discrete allowed states. 
Under Metropolis Monte Carlo dynamics, with random 
initial conditions, we find that at low temperatures the 
model exhibits power law relaxation for the initial stages 
of decay, and at the later stages the relaxation obeys 
a stretched exponential with the exponent (3 ~ 1/4. 
This type of relaxation behaviour is of the Kohlrausch- 
Williams- Watts type as observed experimentally for real 
proteins JT| T |lO| p~2[] . We find that at zero temperature 
the probability distribution function of the energy steps 
encountered along a relaxation path in phase space also 
obeys a stretched exponential form, with another expo- 
nent a ~ 0.39. We show that (3 = a/ (a + 1), yielding 
a value for (3 which is in very good agreement with our 
simulation results. 

We take as our point of departure the model proposed 
by Haliloglu, Bahar, Erman The central idea of this 
model is that all interactions in the protein are governed 
by confining square-law potentials, so all attractions may 
be treated as if the residues interact with each other 
through Hookean forces ||-|| . 

To keep our model very simple, we consider covalent 
bonds as fixed rods of equal length. The residues lo- 
cated at the vertices may be polar P or hydrophobic 
H . All the hydrophobic vertices are to be connected to 
each other with springs of equal stiffness. This feature 
mimicks the effective pressure that is exerted on the hy- 
drophobic residues by the ambient water molecules, and 
results in their being driven to the relatively less exposed 



center of the molecule in the low lying energy states, 
whereas the polar residues are closer to the surface (see 
Fig. 1), a feature that is common to the native configura- 
tions. The constraints placed on the conformations due 
to the rigid chemical bond lengths and restriction of the 
chemical and dihedral angles to discrete values prevent 
the molecule from collapsing to a point. 





FIG. 1. A chain of N = 48 residues, half of which are ran- 
domly chosen to be hydrophobic, (darker beads) shown a) in 
a random initial configuration and b) in a folded state reached 
under Metropolis dynamics. The chain has folded in such a 
way as to leave the polar residues on the outside. (Generated 
using RasMol V2.6) 

It is known that real proteins are distinguished by H- 
P sequences that lead to unique ground states while a 
randomly chosen H-P sequence will typically give rise to 
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a highly degenerate ground state. Nevertheless, in our 
Monte Carlo study we considered a generic H-P sequence 
obtained by choosing fifty percent of the residues to be 
hydrophobic and distributing them randomly along the 
chain. In the absence of detailed knowledge regarding 
the rules singling out the realistic H-P sequences we be- 
lieved this to be in keeping with our statistical approach. 
It might be speculated that the choice of equal probabil- 
ities for encountering H and P groups along the chain, 
and distributing them randomly, maximizes the configu- 
rational entropy of the chain |l3| ] and enhances the "des- 
ignability" giving rise to rather realistic results. 
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FIG. 2. The decay of the energy of an N — 100 chain 
from a random initial configuration, along a zero tempera- 
ture Metropolis trajectory, of 10,000 steps, averaged over 20 
runs. The later stages fit on a stretched exponential curve 
e(t) ~ exp(-t p ) with (3 = 0.234 ± 0.003. The initial stage 
(inset) is fit by a power law e(t) ~ t~ a with a — 0.49 ± 0.01. 

Our model for the protein chain consists of N 
"residues" which are treated as point vertices, connected 
to each other by rigid rods. The "bond angle" q;, at the 
z'th vertex, i = 1, . . . , N — 1, is fixed to be (— l) l a, with 
a = 68°. The dihedral angles fa can take on the values 
of and ±2tt/3. The state (conformation) of the system 
is uniquely specified once the numbers {fa} are given. 
Thus, the residues effectively reside on the vertices of a 
tetrahedral lattice. 

The energy of the molecule is 

i,j i,3 

If we define Qi = 1 for the V th vertex being occupied 
by a hydrophobic residue, and Qi = otherwise, we may 
write Cij — QiQj and 

Vij = [(N h - l)ci,j - Cjj-i - Ci,j+i]Si,j 

- (1 - ^.iX 1 - ■ (2) 

The position vectors r, of each of the vertices in the 
chain can be expressed in terms of a sum over the direc- 
tors Rj of unit length representing the chemical bonds, 



which may be obtained from Ri by successive rotations 
Mfc(afc) and Tk(<j>k) through the bond and the dihedral 
angles @, 

i-l 2 

= Y / Y[ T k(MM k {a k )Ii 1 , (3) 

j=i k=j 

where we may choose Ri along any Cartesian direction 
in our laboratory frame without loss of generality. 
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FIG. 3. The decay of the energy of an N = 48 chain, along 
a Metropolis trajectory, from a random initial configuration 
averaged over 100 runs for 7 = 0.3. The initial stage (inset) 
is fit by a power law e(t) ~ t~ a with a — 0.57 ± 0.01, and the 
late stage to a stretched exponential with /3 = 0.25 ± 0.03. 

In order to investigate the relaxation properties of 
the present model, we have employed Metropolis Monte 
Carlo dynamics. This consisted of a) choosing a pair 
of dihedral angles randomly on the chain, and up- 
dating the (fa , fai ) in a way that preserves angular mo- 
mentum, incrementing them in opposite directions by 
A<p = ±27r/3, b) accepting the move with unit proba- 
bility if AE < and with probability p = exp(— jAE)) 
for AE > 0, c) repeating the second step once before 
discarding the pair altogether and going to the first step. 
Here 7 serves as an effective inverse temperature. Wc 
monitor the relaxation of the total energy as a function 
of "time" measured in the number of MC steps, ( i.e., 
the number of pairs sampled) until a steady state 

is reached, typically in about 10,000 steps. The results 
for chains of N — 100 averaged over 20 randomly chosen 
initial configurations at zero temperature (7 = 00) are 
shown in Fig. 2. Defining e = (E — Eq )/Ei, where Eq is 
the (time- averaged) equilibrium energy and Ej, the ini- 
tial value, we find that it obeys a power law, e(t) ~ t~ a 
with a = 0.49 ± 0.01 for the initial stages of the decay, 
while later stages can be fitted by a stretched exponen- 
tial e(t) - e~ tf3 with (3 = 0.234 ± 0.003. We also per- 
formed simulations for different values of 7, for chains 
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of N = 48, averaging over 100 runs with random initial 
configurations. For 7 — > 00, 7 = 0.5 and 7 = 0.3, the 
above relaxation behaviour continues to hold and the ex- 
ponents do not seem to depend on 7, with (3 ~ 1/4 and 
a ~ 1/2 as given in Table I. 

Table I The exponent a and f3 found for the power law 
and stretched exponential decay of the total energy with time, 
for different chain lengths N and inverse temperatures 7. The 
fits were obtained from a weighted least-squares computation. 
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The variation of the total energy in time is sketched in 
Fig. 4 over a short sequence of relaxation events. Clearly 
one may write E(t), averaged over many independent 
runs, as (E(t)) = (E(0) - J2Zi A E&(t - U)) where 9 
is the Heavyside step function and t 
the time derivative one gets, 
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(E(t)) 



■£Ai^(i-5> fe )) 



(4) 
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At zero temperature, the expectation value of E(t) can 
be calculated by carrying out an integration over the dis- 
tibution of waiting times {Vfe}, and the distribution of 
energy steps encountered along the relaxation path. The 
expectation value, (E(t)) is then, 



A I 



(^))=-(^A^5(t-^T fc )) AE ,r 



(5) 
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FIG. 4. A schematic plot of the variation of the total 
energy with time. 

The distribution of waiting times Tfe is dependent only 
on the configuration of the chain at the K th step and 
independent of the previous waiting times. Since the dy- 
namics is just changing a pair of dihedral angles in oppo- 
site directions, for each conformation {4>i} one may define 



an associated chain of N(N — l)/2 sites, with each site 
corresponding to a pair (i, i') on the original chain. On 
the associated chain, a site will be assigned the value 1 if 
the corresponding pair has at least one "allowed" move, 
and the value if both moves are "blocked." Now the 
probabilities of encountering allowed or blocked moves as 
one implements the Metropolis dynamics outlined above 
are simply given by the density of 1' s or 0' s on the as- 
sociated chain at a given relaxation step, namely, p k and 
qk = 1 — pk- Therefore, in the fc'th conformation, the 
probability of making a transition after Tfe blocked moves 
simply obeys the first passage time distribution |TJ] , 



Pk(Tk) = HkZ 



/i fe = |mg fe | 



(6) 



Writing the 8 function in equation (|^) in the Fourier rep- 
resentation and performing the r-integrals we get 
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We may argue that the larger the energy loss in a relax- 
ation event, the longer it will take for the phase point to 
make a transition out of this state. Since j\ik is roughly 
the expectation for t^, we assume that /ife ~ 1/AEk- 
With the assumption that the energy steps encountered 
along a relaxation path are independently distributed, 
i.e., P(AE 1 . . . AE M ) = nfii P(^Es) for a process of 
M steps, one finds, 



-, m j-i 



(8) 



where Ij/(t) is 



d{AE e ) e &E eP(AE e ) 



n 

fe=0 



AE, 



AEz-AEu, AEk 



(9) 



We have determined from our simulations that the dis- 
tribution P(AEi) (see Fig. 5) also has the stretched 
exponential form P{AEi) = P e~^ El ^ a . The angular 
brackets then take the form 

/>OG 

AE t / (AE e - AE^- 1 exp{-(AE k ) a )dAE k (10) 
Jo 

which we approximate by AE, exp(—(AE k ) a ). The inte- 
gration in equation (|9J) is then straightforward, leading, 
upon substitution in (H), to 
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exp(— Oji' 3 ) 



(11) 
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where a,j = j(l — a)(aj) 13 (1 + 0) 1 and 



a + 1 



(12) 



Substituting the value of a we find from our simulations, 
namely a = 0.39 ± 0.02, we get (3 = 0.28 ± 0.01 which is 
the result we obtained from the fits to the MC simula- 
tions within our error bounds. 
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FIG. 5. The distribution of energy differences encountered 
along the relaxation path are fit to a stretched exponential. 
Level spacing histograms were formed for chains of N = 48 
and averaged over 100 runs for the zero-temperature Metropo- 
lis relaxation. The exponent a of the stretched exponential is 
found to be 0.39 ± 0.02. 

A study of the correlations between fluctuations about 
the native state |l6| for the Beads-and-Springs model, 
with the H-P sequence and the contact map for the native 
states for seven real proteins (61yz, lcd8, lbet, lfil, lbab, 
lcsq and lhiv) was performed by Erman |r^| . Using the 
Gaussian approximation jl8| to the coherent scattering 
function and a normal mode analysis pM E9|j20| , he also 
finds a stretched exponential relaxation with (3 = 1 /4 . 
Experiments on real proteins and polymers [[T 10- 1^] 
yield 0.2 < (3 < 0.4. Our results seem to be closer to 
1/4 and smaller than the values most commonly found 
for spin-glasses namely 1/3. It should also be noted 
that glassy behaviour is obtained here in the absence 
of quenched randomness, or of frustration arising from 
steric hindrances, which we do not take into account. 

Comparing the relaxation behaviour near the native 
state with the behaviour we observe at relatively high 
energies for random heteropolymers, we conclude that 
the relaxation behaviour, and therefore the dynamics and 



the structure of the energy landscape are universal over 
a very large range of energies, and are relatively indepen- 
dent of the specific sequence or the details of the dynam- 
ics. 
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